# load required R libraries
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(ggplot2)
library(AICcmodavg)
# load in data file
# open this file in excel to examine it's structure if you're not familiar with R
pronghorn<-read.csv(file.path("..","pronghorn.csv"), header=TRUE, as.is=TRUE, strip.white=TRUE)
# slope covariate in degrees. Convert to proportion slope
pronghorn$slope<-pronghorn$slope/90
# distance to water in m. Convert to km
pronghorn$DW<-pronghorn$DW/1000
#Set aspect as a factor
pronghorn$aspect = as.factor(pronghorn$aspect)
input.history <- data.frame(freq=1,
ch=apply(pronghorn[,c("Survey.1","Survey.2")],1,paste, collapse=""),
stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 00
## 2 1 01
## 3 1 00
## 4 1 00
## 5 1 01
## 6 1 11
input.history = cbind(input.history, pronghorn[,4:7])
head(input.history)
## freq ch sagebrush slope DW aspect
## 1 1 00 9.0 0.00000000 0.025 W
## 2 1 01 18.0 0.05555556 0.150 S
## 3 1 00 8.4 0.50000000 0.150 W
## 4 1 00 3.2 0.72222222 0.375 E
## 5 1 01 12.0 0.05555556 0.375 S
## 6 1 11 7.8 0.05555556 0.150 S
# Change any NA to . in the chapter history
input.history$ch <- gsub("NA",".", input.history$ch, fixed=TRUE)
head(input.history)
## freq ch sagebrush slope DW aspect
## 1 1 00 9.0 0.00000000 0.025 W
## 2 1 01 18.0 0.05555556 0.150 S
## 3 1 00 8.4 0.50000000 0.150 W
## 4 1 00 3.2 0.72222222 0.375 E
## 5 1 01 12.0 0.05555556 0.375 S
## 6 1 11 7.8 0.05555556 0.150 S
#Create the data structure
prong <- process.data(data = input.history, group = "aspect",
model = "Occupancy")
summary(prong)
## Length Class Mode
## data 7 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 4 data.frame list
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 2 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 4 -none- numeric
## group.covariates 1 data.frame list
## nstrata 1 -none- numeric
## strata.labels 0 -none- NULL
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
###############################
# prepare some data frames that will be used later for creating graphs
sagebrush<-expand.grid(sagebrush=seq(0,20,length.out=100),slope=median(pronghorn$slope),DW=median(pronghorn$DW),aspect=c("N","E","S","W"))
slope<-expand.grid(sagebrush=median(pronghorn$sagebrush),slope=seq(0,1,length.out=100),DW=median(pronghorn$DW),aspect=c("N","E","S","W"))
DW<-expand.grid(sagebrush=median(pronghorn$sagebrush),slope=median(pronghorn$slope),DW=seq(0,3,length.out=100),aspect=c("N","E","S","W"))
aspect<-expand.grid(sagebrush=median(pronghorn$sagebrush),slope=median(pronghorn$slope),DW=median(pronghorn$DW),aspect=c("N","E","S","W"))
#################################################################################################
### logistic regression using GLM
### 'naive' analysis that doesn't account for detection probability
#################################################################################################
glm.data = pronghorn[,4:7]
# y =1 if ever detected at plot, =0 if never detected
glm.data$y<-as.numeric(rowSums(pronghorn[,2:3])>0)
glm.mod1<-glm(y~1,data=glm.data,family=binomial)
glm.mod2<-glm(y~sagebrush,data=glm.data,family=binomial)
glm.mod3<-glm(y~slope,data=glm.data,family=binomial)
glm.mod4<-glm(y~sagebrush+slope,data=glm.data,family=binomial)
glm.mod5<-glm(y~DW,data=glm.data,family=binomial)
glm.mod6<-glm(y~sagebrush+DW,data=glm.data,family=binomial)
glm.mod7<-glm(y~slope+DW,data=glm.data,family=binomial)
glm.mod8<-glm(y~sagebrush+slope+DW,data=glm.data,family=binomial)
glm.mod9<-glm(y~aspect,data=glm.data,family=binomial)
glm.mod10<-glm(y~sagebrush+aspect,data=glm.data,family=binomial)
glm.mod11<-glm(y~slope+aspect,data=glm.data,family=binomial)
glm.mod12<-glm(y~sagebrush+slope+aspect,data=glm.data,family=binomial)
glm.mod13<-glm(y~DW+aspect,data=glm.data,family=binomial)
glm.mod14<-glm(y~sagebrush+DW+aspect,data=glm.data,family=binomial)
glm.mod15<-glm(y~slope+DW+aspect,data=glm.data,family=binomial)
glm.mod16<-glm(y~sagebrush+slope+DW+aspect,data=glm.data,family=binomial)
mods<-list(glm.mod1,glm.mod2,glm.mod3,glm.mod4,
glm.mod5,glm.mod6,glm.mod7,glm.mod8,
glm.mod9,glm.mod10,glm.mod11,glm.mod12,
glm.mod13,glm.mod14,glm.mod15,glm.mod16)
## AIC comparison of GLM models
aictab(mods,second.ord=FALSE)
## Warning in aictab.AICglm.lm(mods, second.ord = FALSE):
## Model names have been supplied automatically in the table
##
## Model selection based on AIC:
##
## K AIC Delta_AIC AICWt Cum.Wt LL
## Mod7 3 351.26 0.00 0.23 0.23 -172.63
## Mod5 2 351.48 0.22 0.21 0.44 -173.74
## Mod8 4 352.08 0.82 0.16 0.60 -172.04
## Mod6 3 352.44 1.18 0.13 0.73 -173.22
## Mod15 6 354.05 2.79 0.06 0.79 -171.03
## Mod13 5 354.34 3.08 0.05 0.84 -172.17
## Mod3 2 355.07 3.81 0.03 0.87 -175.54
## Mod16 7 355.31 4.05 0.03 0.91 -170.65
## Mod14 6 355.71 4.45 0.03 0.93 -171.85
## Mod4 3 355.93 4.67 0.02 0.95 -174.97
## Mod1 1 356.89 5.63 0.01 0.97 -177.45
## Mod11 5 357.37 6.11 0.01 0.98 -173.69
## Mod2 2 357.91 6.65 0.01 0.99 -176.95
## Mod12 6 358.71 7.45 0.01 0.99 -173.36
## Mod9 4 358.93 7.67 0.01 1.00 -175.47
## Mod10 5 360.39 9.13 0.00 1.00 -175.19
summary(mods[[7]]) # top aic-ranked model
##
## Call:
## glm(formula = y ~ slope + DW, family = binomial, data = glm.data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.44059 -1.12575 0.08584 1.11069 1.70310
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6088 0.2369 2.570 0.0102 *
## slope -0.8988 0.6097 -1.474 0.1405
## DW -0.3429 0.1437 -2.387 0.0170 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 354.89 on 255 degrees of freedom
## Residual deviance: 345.26 on 253 degrees of freedom
## AIC: 351.26
##
## Number of Fisher Scoring iterations: 4
#####################################################
# produce some plots of model-averaged estimates to
# illustrate estimated effects
#####################################################
# predict presence probability for different sagebrush values, while fixing other values
glm.ma.sagebrush <- modavgPred(mods,type="response",newdata=sagebrush,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = sagebrush, :
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.sagebrush$mod.avg.pred)
l.se<-glm.ma.sagebrush$uncond.se/(glm.ma.sagebrush$mod.avg.pred*(1-glm.ma.sagebrush$mod.avg.pred))
glm.ma.sagebrush$lower<-plogis(logit-1.96*l.se)
glm.ma.sagebrush$upper<-plogis(logit+1.96*l.se)
#jpeg("PronghornNaiveSagebrush.jpg",width=800,height=800,res=144)
#Convert model averaged results to a data frame
glm.ma.sagebrush = as.data.frame(glm.ma.sagebrush)
#Merge model averaged results with sagebrush data set so that plot can be generated
glm.ma.sagebrush = cbind(glm.ma.sagebrush, sagebrush)
# plot only the esimates for northern aspect
glm.ma.sagebrush = glm.ma.sagebrush[glm.ma.sagebrush$aspect == "N",]
ggplot(data = glm.ma.sagebrush, aes(x = sagebrush, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Sagebrush Density")+
geom_line()+
geom_ribbon(aes(ymin = lower.CL, ymax=upper.CL), alpha = .2)+
xlab("Sagebrush density")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict presence probability for different slope values, while fixing other values
glm.ma.slope <- modavgPred(mods,type="response",newdata=slope,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = slope, :
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.slope$mod.avg.pred)
l.se<-glm.ma.slope$uncond.se/(glm.ma.slope$mod.avg.pred*(1-glm.ma.slope$mod.avg.pred))
glm.ma.slope$lower<-plogis(logit-1.96*l.se)
glm.ma.slope$upper<-plogis(logit+1.96*l.se)
#Convert model averaged results to a data frame
glm.ma.slope = as.data.frame(glm.ma.slope)
#Merge model averaged results with slope data set so that plot can be generated
glm.ma.slope = cbind(glm.ma.slope, slope)
# plot only the esimates for northern aspect
glm.ma.slope = glm.ma.slope[glm.ma.slope$aspect == "N",]
#jpeg("PronghornNaiveSlope.jpg",width=800,height=800,res=144)
ggplot(data = glm.ma.slope, aes(x = slope, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Percent Slope")+
geom_line()+
geom_ribbon(aes(ymin = lower.CL, ymax=upper.CL), alpha = .2)+
xlab("Percent Slope")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict presence probability for different sagebrush values, while fixing other values
glm.ma.DW <- modavgPred(mods,type="response",newdata=DW,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = DW, uncond.se = "revised"):
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.DW$mod.avg.pred)
l.se<-glm.ma.DW$uncond.se/(glm.ma.DW$mod.avg.pred*(1-glm.ma.DW$mod.avg.pred))
glm.ma.DW$lower<-plogis(logit-1.96*l.se)
glm.ma.DW$upper<-plogis(logit+1.96*l.se)
#Convert model averaged results to a data frame
glm.ma.DW = as.data.frame(glm.ma.DW)
#Merge model averaged results with distance from water data set so that plot can be generated
glm.ma.DW = cbind(glm.ma.DW, DW)
# plot only the esimates for northern aspect
glm.ma.DW = glm.ma.DW[glm.ma.DW$aspect == "N",]
#jpeg("PronghornNaiveDW.jpg",width=800,height=800,res=144)
ggplot(data = glm.ma.DW, aes(x = DW, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Distance from Water")+
geom_line()+
geom_ribbon(aes(ymin = lower.CL, ymax=upper.CL), alpha = .2)+
xlab("Distance from Water")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict presence probability for different sagebrush values, while fixing other values
glm.ma.aspect <- modavgPred(mods,type="response",newdata=aspect,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = aspect, :
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.aspect$mod.avg.pred)
l.se<-glm.ma.aspect$uncond.se/(glm.ma.aspect$mod.avg.pred*(1-glm.ma.aspect$mod.avg.pred))
glm.ma.aspect$lower<-plogis(logit-1.96*l.se)
glm.ma.aspect$upper<-plogis(logit+1.96*l.se)
#Convert model averaged results to a data frame
glm.ma.aspect = as.data.frame(glm.ma.aspect)
#Merge model averaged results with distance from water data set so that plot can be generated
glm.ma.aspect = cbind(glm.ma.aspect, aspect)
#jpeg("PronghornNaiveAspect.jpg",width=800,height=800,res=144)
ggplot(data = glm.ma.aspect, aes(x = aspect, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Aspect")+
geom_point()+
geom_errorbar(aes(ymin = lower.CL, ymax=upper.CL), width=.1)+
xlab("Aspect")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))+
scale_x_discrete(labels=c("North","East","South","West"))
#dev.off()
#######################################################
### use occupancy models to fit 16 models for occupancy component
### with general model for detection component
#######################################################
psi.mod1<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~1)
##
## Npar : 8
## -2lnL: 618.2046
## AICc : 634.7876
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.0992008 0.4534794 -0.7896189 0.9880205
## p:sagebrush -0.0083828 0.0306420 -0.0684412 0.0516755
## p:slope -0.5019468 0.5969542 -1.6719771 0.6680834
## p:DW -0.3801192 0.1341593 -0.6430714 -0.1171669
## p:aspectN 1.0302164 0.5337578 -0.0159490 2.0763818
## p:aspectS 0.2186547 0.4556807 -0.6744794 1.1117889
## p:aspectW 0.0325442 0.3870915 -0.7261551 0.7912436
## Psi:(Intercept) 1.2275901 0.4236094 0.3973156 2.0578646
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3739132 0.3739132
## Group:aspectN 0.6259195 0.6259195
## Group:aspectS 0.4263373 0.4263373
## Group:aspectW 0.3815625 0.3815625
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7733965
## Group:aspectN 0.7733965
## Group:aspectS 0.7733965
## Group:aspectW 0.7733965
psi.mod2<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush)
##
## Npar : 9
## -2lnL: 616.6
## AICc : 635.3317
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2327892 0.4651546 -0.6789139 1.1444922
## p:sagebrush -0.0305973 0.0349682 -0.0991350 0.0379403
## p:slope -0.5250128 0.6057552 -1.7122930 0.6622675
## p:DW -0.3939351 0.1345470 -0.6576471 -0.1302230
## p:aspectN 1.0276374 0.5313662 -0.0138404 2.0691152
## p:aspectS 0.2516575 0.4615603 -0.6530007 1.1563158
## p:aspectW 0.0395960 0.3927755 -0.7302440 0.8094359
## Psi:(Intercept) 0.8047650 0.4462285 -0.0698429 1.6793729
## Psi:sagebrush 0.0967337 0.0848505 -0.0695732 0.2630407
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3783969 0.3783969
## Group:aspectN 0.6297815 0.6297815
## Group:aspectS 0.4391283 0.4391283
## Group:aspectW 0.3877542 0.3877542
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7698353
## Group:aspectN 0.7698353
## Group:aspectS 0.7698353
## Group:aspectW 0.7698353
psi.mod3<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~slope),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~slope)
##
## Npar : 9
## -2lnL: 615.4814
## AICc : 634.2131
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.0107202 0.4627442 -0.9176989 0.8962585
## p:sagebrush -0.0089672 0.0306577 -0.0690563 0.0511219
## p:slope 0.7597022 0.9159972 -1.0356524 2.5550569
## p:DW -0.4183130 0.1387979 -0.6903570 -0.1462691
## p:aspectN 0.9547212 0.5427315 -0.1090325 2.0184750
## p:aspectS 0.2711839 0.4769413 -0.6636211 1.2059889
## p:aspectW -0.0070798 0.4001361 -0.7913465 0.7771869
## Psi:(Intercept) 1.6617024 0.5614283 0.5613029 2.7621019
## Psi:slope -2.4427776 1.2304245 -4.8544097 -0.0311456
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3946484 0.3946484
## Group:aspectN 0.6287614 0.6287614
## Group:aspectS 0.4609223 0.4609223
## Group:aspectW 0.3929583 0.3929583
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7651072
## Group:aspectN 0.7651072
## Group:aspectS 0.7651072
## Group:aspectW 0.7651072
psi.mod4<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope)
##
## Npar : 10
## -2lnL: 614.3291
## AICc : 635.2271
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1261223 0.4823440 -0.8192720 1.0715167
## p:sagebrush -0.0277485 0.0353179 -0.0969717 0.0414746
## p:slope 0.6174661 0.9272695 -1.1999822 2.4349144
## p:DW -0.4261022 0.1390519 -0.6986440 -0.1535605
## p:aspectN 0.9714743 0.5461287 -0.0989380 2.0418866
## p:aspectS 0.2913566 0.4824198 -0.6541863 1.2368995
## p:aspectW -0.0024913 0.4059794 -0.7982109 0.7932283
## Psi:(Intercept) 1.2190499 0.5879985 0.0665728 2.3715271
## Psi:sagebrush 0.0786067 0.0770442 -0.0723999 0.2296132
## Psi:slope -2.1743636 1.1973963 -4.5212604 0.1725332
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3996331 0.3996331
## Group:aspectN 0.6374886 0.6374886
## Group:aspectS 0.4711228 0.4711228
## Group:aspectW 0.3990355 0.3990355
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7536614
## Group:aspectN 0.7536614
## Group:aspectS 0.7536614
## Group:aspectW 0.7536614
psi.mod5<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~DW),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~DW)
##
## Npar : 9
## -2lnL: 617.7309
## AICc : 636.4626
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.0412280 0.4701329 -0.8802325 0.9626886
## p:sagebrush -0.0066283 0.0310326 -0.0674522 0.0541956
## p:slope -0.5376998 0.6037819 -1.7211122 0.6457127
## p:DW -0.2823184 0.1988454 -0.6720554 0.1074186
## p:aspectN 1.0401682 0.5485397 -0.0349696 2.1153061
## p:aspectS 0.2127124 0.4699171 -0.7083250 1.1337499
## p:aspectW -0.0115443 0.4056807 -0.8066786 0.7835899
## Psi:(Intercept) 1.4946180 0.6171792 0.2849467 2.7042893
## Psi:DW -0.2878648 0.3913352 -1.0548818 0.4791523
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3894981 0.3894981
## Group:aspectN 0.6435363 0.6435363
## Group:aspectS 0.4410973 0.4410973
## Group:aspectW 0.3867565 0.3867565
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7559078
## Group:aspectN 0.7559078
## Group:aspectS 0.7559078
## Group:aspectW 0.7559078
psi.mod6<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+DW),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + DW)
##
## Npar : 10
## -2lnL: 616.4706
## AICc : 637.3685
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2065638 0.4790686 -0.7324106 1.1455382
## p:sagebrush -0.0291282 0.0359774 -0.0996439 0.0413876
## p:slope -0.5376369 0.6111519 -1.7354947 0.6602209
## p:DW -0.3444901 0.1982931 -0.7331446 0.0441644
## p:aspectN 1.0321802 0.5418322 -0.0298109 2.0941712
## p:aspectS 0.2430861 0.4705463 -0.6791847 1.1653569
## p:aspectW 0.0142206 0.4070578 -0.7836127 0.8120539
## Psi:(Intercept) 0.9625415 0.6498678 -0.3111995 2.2362824
## Psi:sagebrush 0.0884006 0.0861860 -0.0805239 0.2573252
## Psi:DW -0.1540460 0.4146716 -0.9668023 0.6587103
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3878428 0.3878428
## Group:aspectN 0.6400983 0.6400983
## Group:aspectS 0.4468755 0.4468755
## Group:aspectW 0.3912244 0.3912244
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7568626
## Group:aspectN 0.7568626
## Group:aspectS 0.7568626
## Group:aspectW 0.7568626
psi.mod7<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~slope+DW),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
## Warning in file(file, ifelse(append, "a", "w")): cannot open file
## 'C:/Users/Neil/Dropbox/OccupancyCourse-2018-02-12/CourseNotes/
## OccupancySampleData/SingleSpeciesSingleSeason/Pronghorn/RMark
## \markxxx20fc5577783b.tmp': Permission denied
##
## Problem extracting output. Look at MARK output file to see what is wrong.
##
##
## ********Following model failed to run :p(~sagebrush + slope + DW + aspect)Psi(~slope + DW)**********
psi.mod8<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + DW)
##
## Npar : 11
## -2lnL: 614.3291
## AICc : 637.4111
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1261884000 0.4852520 -0.8249054 1.0772823
## p:sagebrush -0.0277517000 0.0354116 -0.0971585 0.0416551
## p:slope 0.6175265000 0.9284583 -1.2022519 2.4373049
## p:DW -0.4262582000 0.1881634 -0.7950585 -0.0574579
## p:aspectN 0.9714659000 0.5461245 -0.0989382 2.0418700
## p:aspectS 0.2913720000 0.4825786 -0.6544822 1.2372262
## p:aspectW -0.0024138000 0.4109612 -0.8078978 0.8030702
## Psi:(Intercept) 1.2186212000 0.6841184 -0.1222509 2.5594932
## Psi:sagebrush 0.0786269000 0.0788198 -0.0758599 0.2331138
## Psi:slope -2.1746803000 1.2247519 -4.5751940 0.2258335
## Psi:DW 0.0005125353 0.4173045 -0.8174043 0.8184293
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3996013 0.3996013
## Group:aspectN 0.6374560 0.6374560
## Group:aspectS 0.4710936 0.4710936
## Group:aspectW 0.3990223 0.3990223
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7537063
## Group:aspectN 0.7537063
## Group:aspectS 0.7537063
## Group:aspectW 0.7537063
psi.mod9<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~aspect)
##
## Npar : 11
## -2lnL: 617.5324
## AICc : 640.6144
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.4003398 0.6482169 -0.8701654 1.6708449
## p:sagebrush -0.0018597 0.0317493 -0.0640883 0.0603689
## p:slope -0.6098451 0.6371491 -1.8586575 0.6389672
## p:DW -0.3800823 0.1375927 -0.6497641 -0.1104006
## p:aspectN 0.7490246 0.8004227 -0.8198040 2.3178532
## p:aspectS -0.3477100 0.8363014 -1.9868608 1.2914407
## p:aspectW -0.2748910 0.6768655 -1.6015473 1.0517654
## Psi:(Intercept) 0.5379061 0.9697562 -1.3628161 2.4386283
## Psi:aspectN 0.6122483 1.2055871 -1.7507025 2.9751990
## Psi:aspectS 1.9200380 4.0922674 -6.1008062 9.9408823
## Psi:aspectW 0.6832722 1.1537158 -1.5780107 2.9445552
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4480957 0.4480957
## Group:aspectN 0.6319652 0.6319652
## Group:aspectS 0.3644555 0.3644555
## Group:aspectW 0.3814829 0.3814829
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6313252
## Group:aspectN 0.7595391
## Group:aspectS 0.9211405
## Group:aspectW 0.7722708
psi.mod10<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + aspect)
##
## Npar : 12
## -2lnL: 616.5109
## AICc : 641.7949
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.3525187 0.6877489 -0.9954691 1.7005065
## p:sagebrush -0.0286137 0.0376700 -0.1024470 0.0452195
## p:slope -0.5125101 0.6278379 -1.7430725 0.7180522
## p:DW -0.3870386 0.1382943 -0.6580955 -0.1159817
## p:aspectN 0.9177164 0.8291509 -0.7074194 2.5428522
## p:aspectS 0.1011425 0.8434129 -1.5519468 1.7542318
## p:aspectW -0.1244842 0.7198525 -1.5353950 1.2864267
## Psi:(Intercept) 0.4703986 1.2185154 -1.9178916 2.8586888
## Psi:sagebrush 0.0923935 0.0927078 -0.0893137 0.2741007
## Psi:aspectN 0.2697689 1.4588602 -2.5895971 3.1291350
## Psi:aspectS 0.3770343 1.6888096 -2.9330325 3.6871011
## Psi:aspectW 0.4228579 1.3816065 -2.2850908 3.1308067
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4116411 0.4116411
## Group:aspectN 0.6365754 0.6365754
## Group:aspectS 0.4363367 0.4363367
## Group:aspectW 0.3818580 0.3818580
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7016046
## Group:aspectN 0.7548631
## Group:aspectS 0.7741658
## Group:aspectW 0.7820765
psi.mod11<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~slope+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~slope + aspect)
##
## Npar : 12
## -2lnL: 615.0065
## AICc : 640.2905
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.0678652 0.6002975 -1.1087180 1.2444484
## p:sagebrush -0.0092784 0.0304535 -0.0689672 0.0504103
## p:slope 0.8581690 0.9086492 -0.9227835 2.6391214
## p:DW -0.4095319 0.1407010 -0.6853058 -0.1337580
## p:aspectN 0.7621498 0.7155689 -0.6403653 2.1646648
## p:aspectS 0.2633683 0.6803600 -1.0701374 1.5968739
## p:aspectW -0.1792780 0.5709985 -1.2984351 0.9398790
## Psi:(Intercept) 1.4152192 1.1464980 -0.8319169 3.6623553
## Psi:slope -2.6247930 1.3941315 -5.3572908 0.1077047
## Psi:aspectN 0.5198959 1.2416696 -1.9137766 2.9535685
## Psi:aspectS 0.0337762 1.1458904 -2.2121690 2.2797214
## Psi:aspectW 0.5221727 1.0941017 -1.6222667 2.6666120
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4206680 0.4206680
## Group:aspectN 0.6087636 0.6087636
## Group:aspectS 0.4858379 0.4858379
## Group:aspectW 0.3777029 0.3777029
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7106570
## Group:aspectN 0.8050989
## Group:aspectS 0.7175525
## Group:aspectW 0.8054559
psi.mod12<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + aspect)
##
## Npar : 13
## -2lnL: 613.8932
## AICc : 641.3973
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1216368 0.5921124 -1.0389034 1.2821771
## p:sagebrush -0.0280938 0.0352284 -0.0971415 0.0409538
## p:slope 0.7438526 0.9090821 -1.0379484 2.5256535
## p:DW -0.4162495 0.1407355 -0.6920911 -0.1404078
## p:aspectN 0.8890925 0.7092774 -0.5010911 2.2792762
## p:aspectS 0.3836901 0.6684401 -0.9264525 1.6938327
## p:aspectW -0.1008680 0.5701466 -1.2183555 1.0166194
## Psi:(Intercept) 1.1729917 1.1174537 -1.0172176 3.3632011
## Psi:sagebrush 0.0817058 0.0811017 -0.0772536 0.2406652
## Psi:slope -2.3882031 1.3081353 -4.9521483 0.1757421
## Psi:aspectN 0.2194644 1.2228207 -2.1772642 2.6161930
## Psi:aspectS -0.2343559 1.1534446 -2.4951074 2.0263957
## Psi:aspectW 0.3069431 1.1035211 -1.8559584 2.4698445
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4071951 0.4071951
## Group:aspectN 0.6256312 0.6256312
## Group:aspectS 0.5020293 0.5020293
## Group:aspectW 0.3830929 0.3830929
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7394258
## Group:aspectN 0.7794477
## Group:aspectS 0.6918176
## Group:aspectW 0.7941181
psi.mod13<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~DW+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~DW + aspect)
##
## Npar : 12
## -2lnL: 616.0833
## AICc : 641.3673
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.3484978000 0.6185548 -0.8638696 1.5608651
## p:sagebrush -0.0007054575 0.0288882 -0.0573262 0.0559153
## p:slope -0.5832516000 0.5785206 -1.7171519 0.5506488
## p:DW -0.1850026000 0.1663117 -0.5109735 0.1409683
## p:aspectN 0.6631636000 0.7592232 -0.8249140 2.1512412
## p:aspectS -0.1603997000 0.7308133 -1.5927938 1.2719944
## p:aspectW -0.7939693000 0.6458990 -2.0599314 0.4719927
## Psi:(Intercept) 1.3514713000 1.0932292 -0.7912580 3.4942006
## Psi:DW -0.7081001000 0.4198861 -1.5310768 0.1148766
## Psi:aspectN 0.5878874000 0.9699113 -1.3131388 2.4889135
## Psi:aspectS 0.5971112000 1.0754675 -1.5108052 2.7050276
## Psi:aspectW 2.6533466000 4.3191414 -5.8121707 11.1188640
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4991674 0.4991674
## Group:aspectN 0.6592229 0.6592229
## Group:aspectS 0.4591587 0.4591587
## Group:aspectW 0.3106034 0.3106034
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6119436
## Group:aspectN 0.7395011
## Group:aspectS 0.7412740
## Group:aspectW 0.9572558
psi.mod14<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+DW+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + DW + aspect)
##
## Npar : 13
## -2lnL: 615.8257
## AICc : 643.3298
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.3558580 0.6322075 -0.8832687 1.5949847
## p:sagebrush -0.0137841 0.0444524 -0.1009108 0.0733425
## p:slope -0.5325788 0.6143066 -1.7366199 0.6714622
## p:DW -0.1946976 0.2066084 -0.5996501 0.2102549
## p:aspectN 0.7445531 0.7776470 -0.7796351 2.2687412
## p:aspectS -0.0339919 0.7732338 -1.5495302 1.4815465
## p:aspectW -0.6026046 0.8425989 -2.2540985 1.0488893
## Psi:(Intercept) 0.9073242 1.3417552 -1.7225160 3.5371643
## Psi:sagebrush 0.0570648 0.1032722 -0.1453487 0.2594784
## Psi:DW -0.5241678 0.6481887 -1.7946176 0.7462821
## Psi:aspectN 0.4672039 0.9786560 -1.4509619 2.3853696
## Psi:aspectS 0.3769413 1.0852734 -1.7501947 2.5040773
## Psi:aspectW 1.4234690 2.4312514 -3.3417838 6.1887218
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4868284 0.4868284
## Group:aspectN 0.6663793 0.6663793
## Group:aspectS 0.4783409 0.4783409
## Group:aspectW 0.3417955 0.3417955
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6181228
## Group:aspectN 0.7208722
## Group:aspectS 0.7023533
## Group:aspectW 0.8704628
psi.mod15<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~slope+DW+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~slope + DW + aspect)
##
## Npar : 13
## -2lnL: 614.3522
## AICc : 641.8563
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1723884 0.6532062 -1.1078957 1.4526725
## p:sagebrush -0.0071742 0.0304443 -0.0668449 0.0524965
## p:slope 0.5381309 1.0974847 -1.6129391 2.6892009
## p:DW -0.2690014 0.2043212 -0.6694709 0.1314680
## p:aspectN 0.6632111 0.7600524 -0.8264916 2.1529138
## p:aspectS 0.1602630 0.7536579 -1.3169064 1.6374325
## p:aspectW -0.5020602 0.7402567 -1.9529634 0.9488430
## Psi:(Intercept) 1.5527438 1.0297020 -0.4654721 3.5709598
## Psi:slope -2.3527025 1.4222704 -5.1403526 0.4349476
## Psi:DW -0.3950014 0.4380304 -1.2535409 0.4635382
## Psi:aspectN 0.6077103 1.0147258 -1.3811524 2.5965730
## Psi:aspectS 0.1325752 0.9687185 -1.7661132 2.0312635
## Psi:aspectW 1.1902794 1.4794450 -1.7094329 4.0899917
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4770346 0.4770346
## Group:aspectN 0.6390596 0.6390596
## Group:aspectS 0.5170775 0.5170775
## Group:aspectW 0.3557207 0.3557207
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6433307
## Group:aspectN 0.7680904
## Group:aspectS 0.6731408
## Group:aspectW 0.8557128
psi.mod16<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 14
## -2lnL: 613.6909
## AICc : 643.4337
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1548598 0.6265924 -1.0732613 1.3829810
## p:sagebrush -0.0235144 0.0371548 -0.0963377 0.0493089
## p:slope 0.6533965 0.9726873 -1.2530706 2.5598637
## p:DW -0.3254091 0.2418093 -0.7993554 0.1485371
## p:aspectN 0.8117553 0.7614303 -0.6806481 2.3041587
## p:aspectS 0.3184463 0.7334922 -1.1191984 1.7560910
## p:aspectW -0.2989206 0.7637047 -1.7957817 1.1979406
## Psi:(Intercept) 1.2758931 1.0928880 -0.8661673 3.4179536
## Psi:sagebrush 0.0679343 0.0851332 -0.0989267 0.2347954
## Psi:slope -2.2922963 1.3402262 -4.9191397 0.3345471
## Psi:DW -0.2490802 0.5168063 -1.2620206 0.7638602
## Psi:aspectN 0.3439041 1.1125470 -1.8366881 2.5244963
## Psi:aspectS -0.1260460 1.0575880 -2.1989186 1.9468266
## Psi:aspectW 0.7050914 1.4224375 -2.0828861 3.4930689
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4437019 0.4437019
## Group:aspectN 0.6423559 0.6423559
## Group:aspectS 0.5230574 0.5230574
## Group:aspectW 0.3716670 0.3716670
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6883202
## Group:aspectN 0.7569779
## Group:aspectS 0.6606589
## Group:aspectW 0.8171818
# compare models by AIC
psi.results<-RMark::collect.models(type = "Occupancy")
# look at AIC table
psi.results
## model
## 10 p(~sagebrush + slope + DW + aspect)Psi(~slope)
## 1 p(~sagebrush + slope + DW + aspect)Psi(~1)
## 11 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope)
## 9 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush)
## 12 p(~sagebrush + slope + DW + aspect)Psi(~DW)
## 13 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + DW)
## 14 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + DW)
## 3 p(~sagebrush + slope + DW + aspect)Psi(~slope + aspect)
## 15 p(~sagebrush + slope + DW + aspect)Psi(~aspect)
## 5 p(~sagebrush + slope + DW + aspect)Psi(~DW + aspect)
## 4 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + aspect)
## 2 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + aspect)
## 7 p(~sagebrush + slope + DW + aspect)Psi(~slope + DW + aspect)
## 6 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + DW + aspect)
## 8 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
## npar AICc DeltaAICc weight Deviance
## 10 9 634.2131 0.0000000 0.258623706 615.4814
## 1 8 634.7876 0.5744286 0.194058530 618.2046
## 11 10 635.2271 1.0139319 0.155774305 614.3291
## 9 9 635.3317 1.1185500 0.147835347 616.6000
## 12 9 636.4626 2.2494600 0.083985497 617.7309
## 13 10 637.3685 3.1553819 0.053393189 616.4706
## 14 11 637.4111 3.1979399 0.052269038 614.3291
## 3 12 640.2905 6.0773333 0.012387743 615.0065
## 15 11 640.6144 6.4012599 0.010535433 617.5324
## 5 12 641.3673 7.1541133 0.007230567 616.0833
## 4 13 641.3973 7.1842049 0.007122592 613.8932
## 2 12 641.7949 7.5817633 0.005838608 616.5109
## 7 13 641.8563 7.6431449 0.005662139 614.3522
## 6 13 643.3298 9.1166849 0.002710225 615.8257
## 8 14 643.4337 9.2205413 0.002573080 613.6909
# regression coefficients for occupancy (psi) and detection (p) from top-ranked model
summary(psi.results[[1]])
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~1)
##
## Npar : 8
## -2lnL: 618.2046
## AICc : 634.7876
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.0992008 0.4534794 -0.7896189 0.9880205
## p:sagebrush -0.0083828 0.0306420 -0.0684412 0.0516755
## p:slope -0.5019468 0.5969542 -1.6719771 0.6680834
## p:DW -0.3801192 0.1341593 -0.6430714 -0.1171669
## p:aspectN 1.0302164 0.5337578 -0.0159490 2.0763818
## p:aspectS 0.2186547 0.4556807 -0.6744794 1.1117889
## p:aspectW 0.0325442 0.3870915 -0.7261551 0.7912436
## Psi:(Intercept) 1.2275901 0.4236094 0.3973156 2.0578646
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.3739132 0.3739132
## Group:aspectN 0.6259195 0.6259195
## Group:aspectS 0.4263373 0.4263373
## Group:aspectW 0.3815625 0.3815625
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7733965
## Group:aspectN 0.7733965
## Group:aspectS 0.7733965
## Group:aspectW 0.7733965
psi.results[[1]]$results$real
## estimate se lcl ucl fixed note
## p gE a0 t1 0.3739132 0.0879090 0.2224447 0.5549139
## p gN a0 t1 0.6259195 0.1013090 0.4174385 0.7962143
## p gS a0 t1 0.4263373 0.0764000 0.2871890 0.5782161
## p gW a0 t1 0.3815625 0.0506148 0.2883666 0.4843770
## Psi gE a0 t1 0.7733965 0.0742394 0.5980425 0.8867399
#####################################################
# produce some plots of model-averaged estimates to
# illustrate estimated effects
#####################################################
# predict occupancy for different sagebrush values with North Aspect,
# while fixing other values
Psi.ma <- RMark::model.average(psi.results, param="Psi")
Psi.ma
## par.index estimate se fixed note group age time
## Psi gE a0 t1 9 0.7587925 0.08983376 E 0 1
## Psi gN a0 t1 10 0.7639931 0.07846864 N 0 1
## Psi gS a0 t1 11 0.7632769 0.08527589 S 0 1
## Psi gW a0 t1 12 0.7670260 0.08114739 W 0 1
## Age Time aspect
## Psi gE a0 t1 0 0 E
## Psi gN a0 t1 0 0 N
## Psi gS a0 t1 0 0 S
## Psi gW a0 t1 0 0 W
ddl = make.design.data(prong)
ddl$Psi # see the index numbers
## par.index model.index group age time Age Time aspect
## 1 1 9 E 0 1 0 0 E
## 2 2 10 N 0 1 0 0 N
## 3 3 11 S 0 1 0 0 S
## 4 4 12 W 0 1 0 0 W
psi.ma.sagebrush <- covariate.predictions(psi.results, indices=10,
data=sagebrush[sagebrush$aspect=="N",])
#jpeg("PronghornPsiSagebrush.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.sagebrush$estimates, aes(x = sagebrush, y = estimate)) +
ggtitle("Model Averaged Occupancy Probability as a function of Sagebrush Density")+
geom_line()+
geom_ribbon(aes(ymin = lcl, ymax=ucl), alpha = .2)+
xlab("Sagebrush density")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict occupancy for different slope values in the North Aspect,
# while fixing other values
psi.ma.slope <- covariate.predictions(psi.results, indices=10,
data=slope[slope$aspect=="N",])
#jpeg("PronghornPsislope.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.slope$estimates, aes(x = slope, y = estimate)) +
ggtitle("Model Averaged Occupancy Probability as a function of Percent Slope")+
geom_line()+
geom_ribbon(aes(ymin = lcl, ymax=ucl), alpha = .2)+
xlab("Percent Slope")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict occupancy for different distance to water values in the North Aspect,
# while fixing other values
psi.ma.DW <- covariate.predictions(psi.results, indices=10,
data=DW[DW$aspect=="N",])
#jpeg("PronghornPsiDW.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.DW$estimates, aes(x = DW, y = estimate)) +
ggtitle("Model Averaged Occupancy Probability as a function of Distance to Water (km)")+
geom_line()+
geom_ribbon(aes(ymin = lcl, ymax=ucl), alpha = .2)+
xlab("Distance to Water (km)")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict occupancy for different aspects, while fixing other values
psi.ma.aspect.E <- covariate.predictions(psi.results, indices=9,
data=aspect[aspect$aspect == "E",])
psi.ma.aspect.N <- covariate.predictions(psi.results, indices=10,
data=aspect[aspect$aspect == "N",])
psi.ma.aspect.S <- covariate.predictions(psi.results, indices=11,
data=aspect[aspect$aspect == "S",])
psi.ma.aspect.W <- covariate.predictions(psi.results, indices=12,
data=aspect[aspect$aspect == "W",])
psi.ma.aspect = rbind(psi.ma.aspect.E$estimates[,7:11],
psi.ma.aspect.N$estimates[,7:11],
psi.ma.aspect.S$estimates[,7:11],
psi.ma.aspect.W$estimates[,7:11])
#jpeg("PronghornPsiaspect.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.aspect, aes(x = aspect, y = estimate)) +
ggtitle("Model Averaged Occupancy Probability as a function of Aspect")+
geom_point()+
geom_errorbar(aes(ymin = lcl, ymax=ucl), width = .1)+
xlab("Aspect")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))+
scale_x_discrete(labels=c("North","East","South","West"))
#dev.off()
##############################################################
### use occupancy models to fit 16 different models for detection component
### with general model for occupancy component
##############################################################
p.mod1<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~1)
)
)
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 8
## -2lnL: 623.3448
## AICc : 639.9278
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1292928 0.1870862 -0.4959817 0.2373961
## Psi:(Intercept) 1.2431420 0.7538943 -0.2344909 2.7207750
## Psi:sagebrush 0.0267133 0.0572365 -0.0854702 0.1388967
## Psi:slope -1.3414061 0.9375843 -3.1790714 0.4962591
## Psi:DW -0.4839058 0.2385330 -0.9514304 -0.0163812
## Psi:aspectN 1.5695425 1.3354381 -1.0479163 4.1870014
## Psi:aspectS 0.3289246 0.7103765 -1.0634134 1.7212626
## Psi:aspectW 0.3410763 0.5886568 -0.8126911 1.4948436
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4677218 0.4677218
## Group:aspectN 0.4677218 0.4677218
## Group:aspectS 0.4677218 0.4677218
## Group:aspectW 0.4677218 0.4677218
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6172738
## Group:aspectN 0.8856984
## Group:aspectS 0.6914523
## Group:aspectW 0.6940387
p.mod2<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 9
## -2lnL: 623.3321
## AICc : 642.0638
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1082644 0.2638714 -0.6254524 0.4089237
## p:sagebrush -0.0040878 0.0362233 -0.0750855 0.0669099
## Psi:(Intercept) 1.2174384 0.7808121 -0.3129533 2.7478302
## Psi:sagebrush 0.0306016 0.0669417 -0.1006042 0.1618074
## Psi:slope -1.3316603 0.9381041 -3.1703444 0.5070237
## Psi:DW -0.4798447 0.2399151 -0.9500784 -0.0096111
## Psi:aspectN 1.5626327 1.3335946 -1.0512127 4.1764780
## Psi:aspectS 0.3235925 0.7096426 -1.0673071 1.7144921
## Psi:aspectW 0.3363162 0.5874769 -0.8151386 1.4877710
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4687213 0.4687213
## Group:aspectN 0.4687213 0.4687213
## Group:aspectS 0.4687213 0.4687213
## Group:aspectW 0.4687213 0.4687213
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6166918
## Group:aspectN 0.8847460
## Group:aspectS 0.6897867
## Group:aspectW 0.6925028
p.mod3<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~slope)
)
)
##
## Output summary for Occupancy model
## Name : p(~slope)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 9
## -2lnL: 623.3045
## AICc : 642.0362
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1567431 0.2301297 -0.6077973 0.2943112
## p:slope 0.1812655 0.8963154 -1.5755127 1.9380438
## Psi:(Intercept) 1.2771137 0.7746926 -0.2412838 2.7955113
## Psi:sagebrush 0.0258068 0.0572145 -0.0863336 0.1379472
## Psi:slope -1.5044923 1.2038988 -3.8641340 0.8551494
## Psi:DW -0.4826341 0.2372009 -0.9475479 -0.0177202
## Psi:aspectN 1.5480278 1.3020080 -1.0039080 4.0999635
## Psi:aspectS 0.3086790 0.7103388 -1.0835851 1.7009431
## Psi:aspectW 0.3442614 0.5868525 -0.8059696 1.4944924
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4697708 0.4697708
## Group:aspectN 0.4697708 0.4697708
## Group:aspectS 0.4697708 0.4697708
## Group:aspectW 0.4697708 0.4697708
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6172047
## Group:aspectN 0.8834720
## Group:aspectS 0.6870534
## Group:aspectW 0.6946526
p.mod4<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+slope)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 10
## -2lnL: 623.2947
## AICc : 644.1927
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1374453 0.3022527 -0.7298605 0.4549699
## p:sagebrush -0.0035748 0.0362319 -0.0745892 0.0674397
## p:slope 0.1749834 0.8995286 -1.5880926 1.9380595
## Psi:(Intercept) 1.2534696 0.8052707 -0.3248610 2.8318002
## Psi:sagebrush 0.0292046 0.0669074 -0.1019340 0.1603431
## Psi:slope -1.4896506 1.2105579 -3.8623442 0.8830429
## Psi:DW -0.4791660 0.2387276 -0.9470721 -0.0112599
## Psi:aspectN 1.5428524 1.3017533 -1.0085842 4.0942890
## Psi:aspectS 0.3047727 0.7096995 -1.0862383 1.6957836
## Psi:aspectW 0.3398694 0.5862654 -0.8092107 1.4889496
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4705636 0.4705636
## Group:aspectN 0.4705636 0.4705636
## Group:aspectS 0.4705636 0.4705636
## Group:aspectW 0.4705636 0.4705636
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6166866
## Group:aspectN 0.8827113
## Group:aspectS 0.6857406
## Group:aspectW 0.6932542
p.mod5<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~DW)
)
)
##
## Output summary for Occupancy model
## Name : p(~DW)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 9
## -2lnL: 618.5398
## AICc : 637.2715
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.3115659 0.2644374 -0.2067314 0.8298632
## p:DW -0.4789734 0.2206550 -0.9114572 -0.0464897
## Psi:(Intercept) 0.8346061 0.8483754 -0.8282097 2.4974218
## Psi:sagebrush 0.0465306 0.0656924 -0.0822266 0.1752877
## Psi:slope -1.7548344 1.1895887 -4.0864283 0.5767596
## Psi:DW 0.2163214 0.6870172 -1.1302325 1.5628752
## Psi:aspectN 2.6868826 5.3632334 -7.8250551 13.1988200
## Psi:aspectS 0.3497866 0.9601012 -1.5320117 2.2315849
## Psi:aspectW 0.0187749 0.8516805 -1.6505189 1.6880687
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4269039 0.4269039
## Group:aspectN 0.4269039 0.4269039
## Group:aspectS 0.4269039 0.4269039
## Group:aspectW 0.4269039 0.4269039
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7224371
## Group:aspectN 0.9745055
## Group:aspectS 0.7869050
## Group:aspectW 0.7261861
p.mod6<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+DW)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + DW)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 10
## -2lnL: 618.3249
## AICc : 639.2228
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.3972993 0.3222992 -0.2344071 1.0290056
## p:sagebrush -0.0159999 0.0346463 -0.0839066 0.0519069
## p:DW -0.4723457 0.2071074 -0.8782763 -0.0664151
## Psi:(Intercept) 0.7385071 0.8371688 -0.9023437 2.3793579
## Psi:sagebrush 0.0623301 0.0745180 -0.0837252 0.2083855
## Psi:slope -1.6620507 1.1119316 -3.8414367 0.5173353
## Psi:DW 0.1749046 0.5744046 -0.9509284 1.3007376
## Psi:aspectN 2.2981358 3.5276458 -4.6160501 9.2123217
## Psi:aspectS 0.3378186 0.9168951 -1.4592958 2.1349330
## Psi:aspectW 0.0421137 0.7896317 -1.5055644 1.5897918
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4336516 0.4336516
## Group:aspectN 0.4336516 0.4336516
## Group:aspectS 0.4336516 0.4336516
## Group:aspectW 0.4336516 0.4336516
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7093234
## Group:aspectN 0.9604652
## Group:aspectS 0.7738044
## Group:aspectW 0.7179294
p.mod7<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~slope+DW)
)
)
##
## Output summary for Occupancy model
## Name : p(~slope + DW)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 10
## -2lnL: 617.9128
## AICc : 638.8107
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2292113 0.2863698 -0.3320734 0.7904961
## p:slope 0.7231122 0.9038455 -1.0484250 2.4946493
## p:DW -0.4923790 0.2029333 -0.8901282 -0.0946298
## Psi:(Intercept) 0.9756900 0.8499409 -0.6901941 2.6415741
## Psi:sagebrush 0.0428947 0.0636151 -0.0817909 0.1675803
## Psi:slope -2.2926527 1.2193173 -4.6825146 0.0972093
## Psi:DW 0.1722928 0.5412005 -0.8884602 1.2330458
## Psi:aspectN 1.9829920 2.2929153 -2.5111221 6.4771062
## Psi:aspectS 0.2373511 0.8841894 -1.4956602 1.9703625
## Psi:aspectW 0.0400283 0.7754901 -1.4799324 1.5599889
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4374592 0.4374592
## Group:aspectN 0.4374592 0.4374592
## Group:aspectS 0.4374592 0.4374592
## Group:aspectW 0.4374592 0.4374592
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7152384
## Group:aspectN 0.9480416
## Group:aspectS 0.7610263
## Group:aspectW 0.7233203
p.mod8<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+slope+DW)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 11
## -2lnL: 617.7068
## AICc : 640.7887
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.3149677 0.3404085 -0.3522330 0.9821684
## p:sagebrush -0.0153588 0.0339322 -0.0818658 0.0511483
## p:slope 0.7105388 0.8958202 -1.0452689 2.4663465
## p:DW -0.4946147 0.1974866 -0.8816884 -0.1075410
## Psi:(Intercept) 0.8824013 0.8528338 -0.7891530 2.5539555
## Psi:sagebrush 0.0574724 0.0718775 -0.0834075 0.1983523
## Psi:slope -2.2166351 1.2078778 -4.5840756 0.1508055
## Psi:DW 0.1658929 0.5011719 -0.8164041 1.1481899
## Psi:aspectN 1.8840624 2.1024749 -2.2367884 6.0049132
## Psi:aspectS 0.2293079 0.8687842 -1.4735092 1.9321249
## Psi:aspectW 0.0410013 0.7529185 -1.4347189 1.5167215
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4415306 0.4415306
## Group:aspectN 0.4415306 0.4415306
## Group:aspectS 0.4415306 0.4415306
## Group:aspectW 0.4415306 0.4415306
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7099649
## Group:aspectN 0.9415456
## Group:aspectS 0.7548264
## Group:aspectW 0.7183344
p.mod9<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 11
## -2lnL: 616.2856
## AICc : 639.3676
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.0692256 0.5554516 -1.0194595 1.1579106
## p:aspectN 0.6446771 0.7241945 -0.7747442 2.0640983
## p:aspectS 0.0836043 0.6790855 -1.2474033 1.4146120
## p:aspectW -0.6313606 0.5955955 -1.7987279 0.5360067
## Psi:(Intercept) 1.5068075 0.8879813 -0.2336358 3.2472509
## Psi:sagebrush 0.0293693 0.0673039 -0.1025464 0.1612849
## Psi:slope -1.6814312 1.0722264 -3.7829950 0.4201327
## Psi:DW -0.6576831 0.2792859 -1.2050834 -0.1102828
## Psi:aspectN 0.5529679 0.8727896 -1.1576998 2.2636356
## Psi:aspectS 0.1355469 0.8200595 -1.4717697 1.7428635
## Psi:aspectW 1.2852118 1.0216623 -0.7172464 3.2876700
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.5172995 0.5172995
## Group:aspectN 0.6712629 0.6712629
## Group:aspectS 0.5381333 0.5381333
## Group:aspectW 0.3630536 0.3630536
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6144127
## Group:aspectN 0.7347519
## Group:aspectS 0.6459871
## Group:aspectW 0.8520929
p.mod10<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 12
## -2lnL: 616.1573
## AICc : 641.4412
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.0938442 0.5561487 -0.9962072 1.1838956
## p:sagebrush -0.0123604 0.0345448 -0.0800681 0.0553473
## p:aspectN 0.6965363 0.7336666 -0.7414501 2.1345228
## p:aspectS 0.1503483 0.6983577 -1.2184329 1.5191294
## p:aspectW -0.5973872 0.6019333 -1.7771765 0.5824021
## Psi:(Intercept) 1.4545937 0.8911207 -0.2920028 3.2011902
## Psi:sagebrush 0.0421738 0.0754768 -0.1057608 0.1901084
## Psi:slope -1.6590106 1.0645884 -3.7456039 0.4275828
## Psi:DW -0.6448484 0.2790734 -1.1918323 -0.0978646
## Psi:aspectN 0.5008700 0.8741650 -1.2124933 2.2142334
## Psi:aspectS 0.0731370 0.8246827 -1.5432411 1.6895150
## Psi:aspectW 1.2216923 1.0239401 -0.7852303 3.2286149
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.5105981 0.5105981
## Group:aspectN 0.6767630 0.6767630
## Group:aspectS 0.5480381 0.5480381
## Group:aspectW 0.3647078 0.3647078
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6195485
## Group:aspectN 0.7287920
## Group:aspectS 0.6366305
## Group:aspectW 0.8467487
p.mod11<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~slope+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~slope + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 12
## -2lnL: 616.2744
## AICc : 641.5584
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.0923294 0.5963840 -1.0765832 1.2612421
## p:slope -0.0938116 0.8850832 -1.8285747 1.6409516
## p:aspectN 0.6436407 0.7250511 -0.7774595 2.0647409
## p:aspectS 0.0655533 0.7015259 -1.3094375 1.4405441
## p:aspectW -0.6474834 0.6153066 -1.8534844 0.5585175
## Psi:(Intercept) 1.4910643 0.9032439 -0.2792938 3.2614224
## Psi:sagebrush 0.0294573 0.0678765 -0.1035807 0.1624953
## Psi:slope -1.5888391 1.4120811 -4.3565181 1.1788399
## Psi:DW -0.6638705 0.2892111 -1.2307242 -0.0970168
## Psi:aspectN 0.5510692 0.8726576 -1.1593399 2.2614782
## Psi:aspectS 0.1548451 0.8470204 -1.5053149 1.8150052
## Psi:aspectW 1.3206065 1.1003370 -0.8360541 3.4772672
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.5184577 0.5184577
## Group:aspectN 0.6720574 0.6720574
## Group:aspectS 0.5347981 0.5347981
## Group:aspectW 0.3604022 0.3604022
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6132321
## Group:aspectN 0.7334090
## Group:aspectS 0.6492544
## Group:aspectW 0.8558851
p.mod12<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+slope+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 13
## -2lnL: 616.1531
## AICc : 643.6572
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1073542 0.5938192 -1.0565314 1.2712399
## p:sagebrush -0.0120981 0.0347851 -0.0802769 0.0560807
## p:slope -0.0571443 0.8807908 -1.7834943 1.6692058
## p:aspectN 0.6948679 0.7347292 -0.7452014 2.1349371
## p:aspectS 0.1382943 0.7240364 -1.2808170 1.5574057
## p:aspectW -0.6074265 0.6224864 -1.8274999 0.6126469
## Psi:(Intercept) 1.4451622 0.9050259 -0.3286886 3.2190130
## Psi:sagebrush 0.0420438 0.0758886 -0.1066978 0.1907854
## Psi:slope -1.6043319 1.3744388 -4.2982319 1.0895682
## Psi:DW -0.6483025 0.2865105 -1.2098632 -0.0867419
## Psi:aspectN 0.5008678 0.8740903 -1.2123493 2.2140849
## Psi:aspectS 0.0852711 0.8500523 -1.5808314 1.7513736
## Psi:aspectW 1.2417834 1.0853784 -0.8855582 3.3691251
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.5114362 0.5114362
## Group:aspectN 0.6771316 0.6771316
## Group:aspectS 0.5458824 0.5458824
## Group:aspectW 0.3631603 0.3631603
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6187042
## Group:aspectN 0.7280833
## Group:aspectS 0.6386069
## Group:aspectW 0.8488789
p.mod13<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~DW + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 12
## -2lnL: 614.3604
## AICc : 639.6444
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2714675 0.6048998 -0.9141362 1.4570711
## p:DW -0.2284710 0.1654887 -0.5528290 0.0958869
## p:aspectN 0.6461819 0.7595197 -0.8424767 2.1348404
## p:aspectS 0.0470496 0.7169251 -1.3581237 1.4522228
## p:aspectW -0.6409565 0.6667167 -1.9477212 0.6658082
## Psi:(Intercept) 1.3355597 1.0564465 -0.7350755 3.4061950
## Psi:sagebrush 0.0373833 0.0744867 -0.1086106 0.1833772
## Psi:slope -1.8403115 1.2032781 -4.1987366 0.5181136
## Psi:DW -0.4628937 0.3994897 -1.2458934 0.3201061
## Psi:aspectN 0.5204664 0.9763326 -1.3931455 2.4340783
## Psi:aspectS 0.1688111 0.9501765 -1.6935350 2.0311571
## Psi:aspectW 1.4657262 1.6372434 -1.7432709 4.6747232
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4955940 0.4955940
## Group:aspectN 0.6521623 0.6521623
## Group:aspectS 0.5073558 0.5073558
## Group:aspectW 0.3410585 0.3410585
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6325582
## Group:aspectN 0.7433920
## Group:aspectS 0.6708470
## Group:aspectW 0.8817319
p.mod14<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 13
## -2lnL: 614.1157
## AICc : 641.6199
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.3014115 0.6067217 -0.8877630 1.4905860
## p:sagebrush -0.0181787 0.0383418 -0.0933286 0.0569712
## p:DW -0.2466428 0.1953964 -0.6296198 0.1363342
## p:aspectN 0.7380382 0.7813021 -0.7933140 2.2693904
## p:aspectS 0.1513885 0.7412327 -1.3014276 1.6042046
## p:aspectW -0.5273634 0.7632886 -2.0234091 0.9686824
## Psi:(Intercept) 1.1852688 1.1317635 -1.0329878 3.4035253
## Psi:sagebrush 0.0605010 0.0894744 -0.1148688 0.2358708
## Psi:slope -1.7474230 1.1934335 -4.0865528 0.5917068
## Psi:DW -0.3914992 0.4875408 -1.3470792 0.5640808
## Psi:aspectN 0.4299441 1.0035752 -1.5370634 2.3969515
## Psi:aspectS 0.0756846 0.9702664 -1.8260376 1.9774069
## Psi:aspectW 1.1541047 1.7035039 -2.1847630 4.4929723
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4784296 0.4784296
## Group:aspectN 0.6573941 0.6573941
## Group:aspectS 0.5162576 0.5162576
## Group:aspectW 0.3512160 0.3512160
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6451427
## Group:aspectN 0.7364697
## Group:aspectS 0.6622731
## Group:aspectW 0.8521873
p.mod15<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~slope + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 13
## -2lnL: 614.0996
## AICc : 641.6037
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1360603 0.6524113 -1.1426659 1.4147865
## p:slope 0.5449957 1.0551077 -1.5230154 2.6130068
## p:DW -0.2832781 0.2178792 -0.7103214 0.1437651
## p:aspectN 0.6862052 0.7558291 -0.7952199 2.1676303
## p:aspectS 0.1682808 0.7345743 -1.2714849 1.6080465
## p:aspectW -0.4561378 0.7471015 -1.9204568 1.0081812
## Psi:(Intercept) 1.3941140 1.0666510 -0.6965220 3.4847499
## Psi:sagebrush 0.0391434 0.0710649 -0.1001439 0.1784307
## Psi:slope -2.3030260 1.4016694 -5.0502980 0.4442460
## Psi:DW -0.3496460 0.4649426 -1.2609336 0.5616416
## Psi:aspectN 0.4887017 1.0524039 -1.5740099 2.5514134
## Psi:aspectS 0.0301771 1.0105080 -1.9504185 2.0107728
## Psi:aspectW 1.0204662 1.4521883 -1.8258229 3.8667553
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4712547 0.4712547
## Group:aspectN 0.6390163 0.6390163
## Group:aspectS 0.5132900 0.5132900
## Group:aspectW 0.3609494 0.3609494
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6595538
## Group:aspectN 0.7595124
## Group:aspectS 0.6662969
## Group:aspectW 0.8431391
p.mod16<-RMark::mark(prong,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~sagebrush+slope+DW+aspect),
p =list(formula=~sagebrush+slope+DW+aspect)
)
)
##
## Output summary for Occupancy model
## Name : p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 14
## -2lnL: 613.6909
## AICc : 643.4337
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.1548592 0.6265928 -1.0732628 1.3829811
## p:sagebrush -0.0235144 0.0371547 -0.0963375 0.0493087
## p:slope 0.6534028 0.9726881 -1.2530659 2.5598715
## p:DW -0.3254097 0.2418094 -0.7993563 0.1485368
## p:aspectN 0.8117538 0.7614305 -0.6806501 2.3041576
## p:aspectS 0.3184460 0.7334930 -1.1192004 1.7560924
## p:aspectW -0.2989204 0.7637067 -1.7957856 1.1979448
## Psi:(Intercept) 1.2758915 1.0928892 -0.8661713 3.4179543
## Psi:sagebrush 0.0679346 0.0851331 -0.0989264 0.2347955
## Psi:slope -2.2923018 1.3402267 -4.9191463 0.3345426
## Psi:DW -0.2490792 0.5168062 -1.2620194 0.7638609
## Psi:aspectN 0.3439081 1.1125508 -1.8366914 2.5245076
## Psi:aspectS -0.1260448 1.0575941 -2.1989292 1.9468397
## Psi:aspectW 0.7050933 1.4224455 -2.0828998 3.4930865
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4437019 0.4437019
## Group:aspectN 0.6423555 0.6423555
## Group:aspectS 0.5230573 0.5230573
## Group:aspectW 0.3716670 0.3716670
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6883201
## Group:aspectN 0.7569786
## Group:aspectS 0.6606590
## Group:aspectW 0.8171820
#Delete old occupancy models before creating new AIC table
rm(psi.mod1,psi.mod2,psi.mod3,psi.mod4,psi.mod5,psi.mod6,psi.mod7,psi.mod8,
psi.mod9,psi.mod10,psi.mod11,psi.mod12,psi.mod13,psi.mod14,psi.mod15,psi.mod16)
# compare models by AIC
p.results<-RMark::collect.models(type = "Occupancy")
# look at AIC table
p.results
## model
## 12 p(~DW)Psi(~sagebrush + slope + DW + aspect)
## 14 p(~slope + DW)Psi(~sagebrush + slope + DW + aspect)
## 13 p(~sagebrush + DW)Psi(~sagebrush + slope + DW + aspect)
## 16 p(~aspect)Psi(~sagebrush + slope + DW + aspect)
## 5 p(~DW + aspect)Psi(~sagebrush + slope + DW + aspect)
## 1 p(~1)Psi(~sagebrush + slope + DW + aspect)
## 15 p(~sagebrush + slope + DW)Psi(~sagebrush + slope + DW + aspect)
## 2 p(~sagebrush + aspect)Psi(~sagebrush + slope + DW + aspect)
## 3 p(~slope + aspect)Psi(~sagebrush + slope + DW + aspect)
## 7 p(~slope + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
## 6 p(~sagebrush + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
## 10 p(~slope)Psi(~sagebrush + slope + DW + aspect)
## 9 p(~sagebrush)Psi(~sagebrush + slope + DW + aspect)
## 8 p(~sagebrush + slope + DW + aspect)Psi(~sagebrush + slope + DW + aspect)
## 4 p(~sagebrush + slope + aspect)Psi(~sagebrush + slope + DW + aspect)
## 11 p(~sagebrush + slope)Psi(~sagebrush + slope + DW + aspect)
## npar AICc DeltaAICc weight Deviance
## 12 9 637.2715 0.000000 0.269907202 618.5398
## 14 10 638.8107 1.539222 0.125019193 617.9128
## 13 10 639.2228 1.951352 0.101738141 618.3249
## 16 11 639.3676 2.096110 0.094634609 616.2856
## 5 12 639.6444 2.372883 0.082404212 614.3604
## 1 8 639.9278 2.656319 0.071515834 623.3448
## 15 11 640.7887 3.517240 0.046500276 617.7068
## 2 12 641.4412 4.169733 0.033555861 616.1573
## 3 12 641.5584 4.286853 0.031647259 616.2744
## 7 13 641.6038 4.332255 0.030936934 614.0996
## 6 13 641.6199 4.348375 0.030688585 614.1157
## 10 9 642.0362 4.764690 0.024921558 623.3045
## 9 9 642.0638 4.792280 0.024580125 623.3321
## 8 14 643.4337 6.162181 0.012391213 613.6909
## 4 13 643.6572 6.385705 0.011080931 616.1531
## 11 10 644.1927 6.921192 0.008478065 623.2947
# regression coefficients for occupancy (psi) and detection (p) from top-ranked model
summary(p.results[[1]])
## Output summary for Occupancy model
## Name : p(~1)Psi(~sagebrush + slope + DW + aspect)
##
## Npar : 8
## -2lnL: 623.3448
## AICc : 639.9278
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1292928 0.1870862 -0.4959817 0.2373961
## Psi:(Intercept) 1.2431420 0.7538943 -0.2344909 2.7207750
## Psi:sagebrush 0.0267133 0.0572365 -0.0854702 0.1388967
## Psi:slope -1.3414061 0.9375843 -3.1790714 0.4962591
## Psi:DW -0.4839058 0.2385330 -0.9514304 -0.0163812
## Psi:aspectN 1.5695425 1.3354381 -1.0479163 4.1870014
## Psi:aspectS 0.3289246 0.7103765 -1.0634134 1.7212626
## Psi:aspectW 0.3410763 0.5886568 -0.8126911 1.4948436
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4677218 0.4677218
## Group:aspectN 0.4677218 0.4677218
## Group:aspectS 0.4677218 0.4677218
## Group:aspectW 0.4677218 0.4677218
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6172738
## Group:aspectN 0.8856984
## Group:aspectS 0.6914523
## Group:aspectW 0.6940387
p.results[[1]]$results$real
## estimate se lcl ucl fixed note
## p gE a0 t1 0.4677218 0.0465766 0.3784854 0.5590719
## Psi gE a0 t1 0.6172738 0.1297775 0.3546426 0.8255886
## Psi gN a0 t1 0.8856984 0.1318783 0.3762017 0.9900558
## Psi gS a0 t1 0.6914523 0.1094641 0.4504819 0.8596708
## Psi gW a0 t1 0.6940387 0.0832344 0.5126992 0.8302402
#####################################################
# produce some plots of model-averaged estimates to
# illustrate estimated effects
#####################################################
# predict detection probability for different sagebrush values with North Aspect,
# while fixing other values
p.ma <- RMark::model.average(p.results, param="p")
p.ma
## par.index estimate se fixed note group age time Age
## p gE a0 t1 1 0.4589946 0.10005787 E 0 1 0
## p gE a1 t2 2 0.4589946 0.10005787 E 1 2 1
## p gN a0 t1 3 0.5118833 0.12964052 N 0 1 0
## p gN a1 t2 4 0.5118833 0.12964052 N 1 2 1
## p gS a0 t1 5 0.4675352 0.08433620 S 0 1 0
## p gS a1 t2 6 0.4675352 0.08433620 S 1 2 1
## p gW a0 t1 7 0.4119135 0.06973346 W 0 1 0
## p gW a1 t2 8 0.4119135 0.06973346 W 1 2 1
## Time aspect
## p gE a0 t1 0 E
## p gE a1 t2 1 E
## p gN a0 t1 0 N
## p gN a1 t2 1 N
## p gS a0 t1 0 S
## p gS a1 t2 1 S
## p gW a0 t1 0 W
## p gW a1 t2 1 W
ddl = make.design.data(prong)
ddl$p # see the index numbers
## par.index model.index group age time Age Time aspect
## 1 1 1 E 0 1 0 0 E
## 2 2 2 E 1 2 1 1 E
## 3 3 3 N 0 1 0 0 N
## 4 4 4 N 1 2 1 1 N
## 5 5 5 S 0 1 0 0 S
## 6 6 6 S 1 2 1 1 S
## 7 7 7 W 0 1 0 0 W
## 8 8 8 W 1 2 1 1 W
p.ma.sagebrush <- covariate.predictions(p.results, indices=3,
data=sagebrush[sagebrush$aspect=="N",])
#jpeg("PronghornPSagebrush.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.sagebrush$estimates, aes(x = sagebrush, y = estimate)) +
ggtitle("Model Averaged Detection Probability as a function of Sagebrush Density")+
geom_line()+
geom_ribbon(aes(ymin = lcl, ymax=ucl), alpha = .2)+
xlab("Sagebrush density")+
ylab("Model averaged detection probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict detection probability for different slope values in the North Aspect,
# while fixing other values
p.ma.slope <- covariate.predictions(p.results, indices=3,
data=slope[slope$aspect=="N",])
#jpeg("PronghornPslope.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.slope$estimates, aes(x = slope, y = estimate)) +
ggtitle("Model Averaged Detection Probability as a function of Percent Slope")+
geom_line()+
geom_ribbon(aes(ymin = lcl, ymax=ucl), alpha = .2)+
xlab("Percent Slope")+
ylab("Model averaged detection probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict detection probability for different distance to water values in the North Aspect,
# while fixing other values
p.ma.DW <- covariate.predictions(p.results, indices=3,
data=DW[DW$aspect=="N",])
#jpeg("PronghornPDW.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.DW$estimates, aes(x = DW, y = estimate)) +
ggtitle("Model Averaged Detection Probability as a function of Distance to Water (km)")+
geom_line()+
geom_ribbon(aes(ymin = lcl, ymax=ucl), alpha = .2)+
xlab("Distance to Water (km)")+
ylab("Model averaged detection probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict detection probabilities for different aspects, while fixing other values
p.ma.aspect.E <- covariate.predictions(p.results, indices=1,
data=aspect[aspect$aspect == "E",])
p.ma.aspect.N <- covariate.predictions(p.results, indices=3,
data=aspect[aspect$aspect == "N",])
p.ma.aspect.S <- covariate.predictions(p.results, indices=5,
data=aspect[aspect$aspect == "S",])
p.ma.aspect.W <- covariate.predictions(p.results, indices=7,
data=aspect[aspect$aspect == "W",])
p.ma.aspect = rbind(p.ma.aspect.E$estimates[,7:11],
p.ma.aspect.N$estimates[,7:11],
p.ma.aspect.S$estimates[,7:11],
p.ma.aspect.W$estimates[,7:11])
#jpeg("PronghornPaspect.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.aspect, aes(x = aspect, y = estimate)) +
ggtitle("Model Detection Occupancy Probability as a function of Aspect")+
geom_point()+
geom_errorbar(aes(ymin = lcl, ymax=ucl), width = .1)+
xlab("Aspect")+
ylab("Model averaged detection probability")+
ylim(c(0,1))+
scale_x_discrete(labels=c("North","East","South","West"))
#dev.off()
#cleanup
cleanup(ask=FALSE)